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ABSTRACT 

Recent developments in instrumentation (e.g., in particular the Kepler and CoRoT satellites) pro- 
vide a new opportunity to improve the models of stellar pulsations. Surface layers, rotation, and 
magnetic fields imprint erratic frequency shifts, trends, and other non-random behavior in the fre- 
quency spectra. As our observational uncertainties become smaller, these are increasingly important 
and difficult to deal with using standard fitting techniques. To improve the models, new ways to com- 
pare their predictions with observations need to be conceived. In this paper we present a completely 
probabilistic (Bayesian) approach to asteroseismic model fitting. It allows for varying degrees of prior 
mode identification, corrections for the discrete nature of the grid, and most importantly implements 
a treatment of systematic errors, such as the "surface effects." It removes the need to apply semi- 
empirical corrections to the observations prior to fitting them to the models and results in a consistent 
set of probabilities with which the model physics can be probed and compared. As an example, we 
show a detailed asteroseismic analysis of the Sun. We find a most probable solar age, including a 
35 ± 5 million year pre-main sequence phase, of 4.591 billion years, and initial element mass fractions 
of X = 0.72, Y = 0.264, Z = 0.016, consistent with recent asteroseismic and non-asteroseismic 
studies. 



1. INTRODUCTION 

The success of recent space missions CoRoT and Ke- 
pler, designed for the discovery of exoplanets and the 
analysis of stellar pulsation, have produced a large num - 
ber of high-quality light curves (Chaplin et al. 2010). 
With these data sets, obtained over long time bases of 
several months, we are able to detect variability with 
semi-amplitudes down to a few parts per million. These 
observations have now firmly established the existence 
of solar-type pulsation in a large number of solar-like 
and red-giant stars. Moreover, observations of an un- 
precedented number of <5Scuti stars and other types of 
pulsators have also revealed rich mode spectra. 

These data are now causing a paradigm shift for many 
topics in stellar astrophysics. In particular, the determi- 
nation of fundamental stellar parameters, and any infer- 
ences regarding the physics of stellar interiors, have for 
a long time been restricted to testing theoretical models 
using classic observables such as photometric indices or 
spectroscopic data. Even though these methods have be- 
come more advanced, for instance by applyi ng complex 
Bayesian methods to determine ste llar ages ( Pont & Eyer 
2004 J0rgensen &: Linde grcn 200 51 and to evaluate com - 
peting models (T'akeda et al. 2007 Bazot et al. 2008), 



the value of additional information provided by pulsa- 
tion modes is tremendous, as they directly probe the 
whole star. Already, the asteroseismic community is suc- 
cessful in extracting general characteristics of t he mode 
spectra for many different types o f stars (e.g 



et al. 



2010; 



Mathur 



Kallinger et al. 2010c) and also in devising 
promising tools fo r a co mparati ve interpretat ion of the 
observations (e.g., Bedding & Kjeldsen 2010). Average 



mode parameters, such as the large and small frequency 
separations, and the frequency of maximum power, have 
been shown to successfully constrain stellar parameters 
although certain co rrelations remain as a source for un- 
certainty (see, e.g., |Kallinger et al. 2010b Huber et al 



2011 



into t 



Gai et al. 2011| ). These have been incorporated 



le current advanced p robabilistic pipeline s to in- 



vestigate stellar model grids ( Quirion et al. | 2010|) and al 



ready been applied to recent observations ( |Metcalfe et al 
2010|). The next step to improving our knowledge about 



stellar interiors is to analyze individual pulsation modes 
in an equally rigorous way, to see where our models agree 
or disagree 



In the past, \ -minimization techniques (Guenther 
|fc Bro wn 2004]), or equ ivalent Bayesian analyses (e.g., 
IKallinger et al. 2010a), have been introduced to find 



the pulsation model that most closely reproduce the 
observed frequencies within a large and dense grid 
of models. The Bayesian analysis, in this context, 
only provides an additional framework for constraining 
solutions to models that match our prior knowledge 
about the stars' fundamental parameters. Due to the 
rich information provided by the pulsation frequencies, 
these approaches should be successful in many cases, 
which is why they are being applied al so to the most 
recent Kepler data sets. For instance, [Metcalfe et al.| 
(20101 test various approaches from different modelers 
with different methods that actually use the individual 
frequencies. However, there are currently (at least) 
three major problems when applying these techniques. 

Stellar rotation, at all but the slowest rotation speeds, 
has been shown to produce rotational splittings which 
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Figure 1. E chelle diagram of solar p modes taken from |Broomhall1 
|et al.| 112009 ) (filled circles) and an appropriate solar model con- 
structed using YREC. The higher order model frequencies are in- 
creasingly deviating from the observations due to deficiencies in 
modeling the upper stellar layers. The systematic errors of the 
models are much bigger than the random observational uncertain- 
ties. 



are incompatible with the traditional linear approxima- 
tions. It even perturbs the va lues of the axisymmetric 



[m = 0) frequencies (e.g., see Deupree & Beslin 2010 
and references therein). In order to correctly take this 
into account, the rotation speed as a function of stellar 
depth needs to be known, and extensive computations 
would be necessary to do these effects justice. Given the 
large variety of possible rotation profile characteristics, 
this would greatly expand the dimensionality and size 
of the pulsation model grid. This implies currently 
insurmountable computational expenses for the types 
and sizes of grids that are necessary for a comprehensive 
asteroseismic analysis of many stars. 

For stars with a convective envelope, model frequencies 
at high radial orders differ from observations due to prob- 
lems in modeling the outer layers (see Figure \U. These 
so-called surface effects can be c ompensated by l ooking 
at ratios of frequency differences ( Roxburgh 2005 ) , or by 
"correcting" the observed frequencies through calibra- 
tio n of the surface effects seen for the Sun as proposed 
by Kjeldsen et al.| ( |2008[ ). It is likely that the surface 
correction as calibrated tor the Sun is not universally ap- 
plicable, and eviden ce for this has been mounting (e.g., 
Bedding et al.|2010| ). Moreover, neglecting (or correcting 
for) the surface effects in the observed frequencies is only 
reasonable when studying properties of the star for which 
the outer layers are unimportant. However, if we want 
the theoretical models to more closely reflect reality, we 
need to include more and better physics to bring the com- 
puted frequencies closer to the (un-corrected) observed 
ones. 



Furthermore, the fact that static asteroseismic grids 
can only have a finite resolution in parameter space is 
often neglected. If the error bars of the observed fre- 
quencies are small compared to the differences between 
calculated frequencies in adjacent grid points, the likeli- 
hood of having a model in the grid that corresponds to 
the best model one's code could deliver decreases rapidly. 
The problem of finding the "true" model and the actual 
uncertainties with respect to the grid becomes apparent. 
Even grids with adaptive resolution have the same prob- 
lem in principle, as the decision for further refining the 
resolution of a particular region in parameter space must 
always depend on a number of discrete grid points. This 
problem is much more severe if our aim is to calculate 
probabilities (or some summary statistics) to compare 
different model grids. 

In this paper we present a new approach to asteroseis- 
mic model grid fitting. Our goal is to find a new way 
of putting our model physics to the test that can handle 
all of the aforementioned difficulties. Even restricted to 
models that are unable to produce all the details of the 
observations, we want to know which models are most 
"correct" (i.e., consistent with appropriate fundamental 
parameters and physics), and how well the solution is 
constrained. We show how to quantitatively assess our 
model grids as a function of the observational uncertain- 
ties, the uncertainties of the calculated frequencies, and 
our general prior knowledge about the star and possible 
shortcomings of our models. 

2. BAYESIAN TREATMENT OF SYSTEMATIC ERRORS 

2.1. Basics of Bayesian inference 

Bayes' theorem, applied to the problem of inference, 
states that the probability of a particular hypothesis af- 
ter obtaining new data (i.e., the posterior) is proportional 
to the probability of the hypothesis prior to obtaining the 
new data (i.e., the prior) times the likelihood of obtaining 
the new data, under the assumption that the hypothe- 
sis is true (i.e., the likelihood function). This approach 
to inference is derived from the product and sum rules 
of probability theory that have shown to be necessary 
and suffic ient for consistent, quant itative logical reason- 
im£Jsee Jaynes fc Bretthorst||2003 ). 

m this paper, we stay as close as pos s ible to th e general 
notati on used in Jaynes & Bretthorst ( 2003 ) or Gregory 
( |2005| . We start with Bayes' theorem applied to the 
problem of comparing observations with the predictions 
of a model M. If the predictions of a model M are gov- 
erned by a set of n parameters = {9i, ...,#„}, and we 
define the observations to be represented by the symbol 
D (for data), it is commonly formulated by expressing 
the posterior probability 



P(0\M,D,I) 



P(0\M,I)P(D\9,M,I) 
P(D\M,I) ' 



(1) 



The symbol / is equivalent to the prior information about 
the problem that is investigated. The first term in the 

1 Strictly speaking, Bayes' theorem is only one result that derives 
from these rules. Consistent use of Bayes' theorem, in particular 
the assignment of the various terms in Equation Q , also requires 
knowledge of its origin and consistent application of the product 
and sum rule. However, for the sake of brevity we will simply 
call our approach in this manuscript to be "Bayesian" rather than 
"based on probability theory as extended logic" . 
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numerator of Equation ([I]) is the prior probability of a 
particular set of parameter values 9, given the model 
M and our prior information / about the problem. It 
is independent of any new data which are supposed to 
be analyzed. The second term in the numerator is called 
the likelihood. It gives the likelihood of obtaining the ob- 
served data under the assumption that the predictions of 
model M are correct, given the particular choice of its 
parameter values 0. The denominator in Equation ([TJ is 
called the global likelihood, or evidence, and is the sum 
(or integral) of the numerator over the whole parameter 
space of model M. It therefore acts as a normalization 
constant. Most importantly, if the prior probabilities are 
adequately normalized, it also represents the likelihood 
of obtaining the data given the whole model M , indepen- 
dent of the particular choice of 0. Thus, it can be used as 
a likelihood for comparisons among different alternative 
models. 

More details on the application of Bayes' theorem, in 
particular with respect to dat a analysis in astronomy, 
can be found in Gregory ( 2005 ) . 



d6„ 



2.2. Systematic errors in the Bayesian framework 

One of the strengths of the Bayesian framework is that 
a parameter 6 n , known to be necessary to describe a 
model M, can be marginalized by applying the sum rule. 
In the case of continuous parameters, the sum turns into 
an integral, and by integrating the full posterior over the 
parameter range of n , one obtains the marginal posterior 

P(6 1 ,...,6 n - 1 \M,D,I) = J P(0 1 ,...,0 n _ 1 ,0 n |M,D,/) 

(2) 

The marginal posterior retains the overall effects of in- 
cluding parameter 9 n in the model, but is independent of 
any particular choice of its value. In other words, 9 n is 
"removed" from the detailed analysis. This is similar to 
what is done for calculating the evidence in the denom- 
inator in Equation ([l]) . The only difference is that the 
evidence is the marginal likelihood over all parameters 
of the model, weighted by the prior. 

The reason this is useful is that if the data and the 
model are known to show systematic differences, like 
shifts or trends, such "systematic errors" can simply 
be encoded introducing additional parameters to the 
model M, so that M is able to model these effects as 
well. By subsequently marginalizing over these "fudge 
parameters" , one is then able to perform a standard 
Bayesian analysis without any need for knowing the 
exact value of the systematic error(s). However, even 
though the exact value is unknown, the presence of 
the error is being considered in the evaluation of the 
posterior probabilities. Furthermore, an increasing 
number of "fudge parameters" comes at a cost, because 
it potentially decreases the evidence for the model 
due to the increase in prior volume. As mentioned in 
Section l2~T1 the evidence is used as the value for the 
likelihood of obtaining the data in Bayesian model 
comparison. It is therefore possible to compare models 
with and without "fudge parameters" . Improved models 
that do not need them, but are able to explain the 
observations just as well, will be favored. 



3. TOWARD A BAYESIAN SOLUTION TO 
ASTEROSEISMIC MODEL FITTING 

3.1. Review and problems of the standard approach 

The general problem of asteroseismic model fitting is 
to match observed frequencies fi :0 to those calculated 
from models fi. m . If the n Q b s observed frequencies have 
individual uncertainties cr.; i0 , and the model frequencies 
have random uncertainties <7j !m then a x 2 ~statistic can 
be calculated according to 



2 1 \ (/i,o fi.m) 



^obs 



i=l 



o~7 



(3) 



Searching a large grid of N stellar models Mj with fun- 
damental parameters close to those estimated for the 
observed star will produce a minimum in \ 2 ( = best- 
fit model). In addition, uncertainties can be estimated 
from the change in % 2 as the distance in parameter space 
to the best fit increases. Calculated with adequate stel- 
lar evolution and pulsation codes, it should be possible 
to infer details about the stellar interior and to obtain 
precise fundamental parameters. 

In order to consistently encode prior information about 
the fundamental parameters and other model properties, 
and to make use of all the additional advantages that 
come with the Bayesian approach (all of which will be- 
come clear in the next section), it is much easier to per- 
form the model fitting using probabilities. Assuming that 
the random uncertainties (7i i0 and a^ m are compatible 
with Normal distributions, one can define 

°1 = °f,o + (4) 

This leads to the likelihood for observing the data (= 
the specific values of /i j0 ), given a single observed and 
calculated frequency 



i,oi— j 



Mj,I) 



1 



27TCT, 



exp 



{fi,o fi,in) 

2^| 



(5) 

Here, /,:, M.j,m stands for the proposition 11 The observed 
mode fi j0 corresponds to the calculated mode /i.m-'j^Nat- 
urally, we want our models Mj to reproduce all observed 
frequencies. Assuming that each observed frequency is a 
statistically independent datapoint, this leads to a prod- 
uct for the likelihood of obtaining all observed frequency 
values given that the model is correct 



I'D Mj,I) = J! P(/ i)0 |/i,»4i,m,Af J , J). 



(6) 



Here, D stands for complete set of observed frequencies 
and their uncertainties. This can then be incorporated 
in the usual framework for Bayesian inference. 

Alas, both of the mentioned, straightforward ap- 
proaches above suffer from the following problems: 

2 Although the explicit notation seems clumsy at first glance, it 
is actually one of the major assets of the Bayesian approach. It 
visualizes exactly which propositions we are evaluating, and under 
which conditions the probabilities are calculated. Slightly different 
propositions or conditions can yield vastly different results. If the 
notation is explicit, there are no hidden variables or assumptions. 
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1. The most appropriate model is not necessarily 
the one that minimizes Equation ^ or maximizes 
Equation There are many possible scenarios 
where this would be the case (e.g., due to surface 
effects, stellar activity, magnetic field effects, rota- 
tional effects). Straightforward application of the 
formalism above will then lead to wrong or non- 
sensical results in both best fit and derived uncer- 
tainties. Even worse, this would propagate into our 
assessment of the model physics that were used to 
produce the models. 

2. In case of such systematic differences, we need to 
take into account that for each observed pulsation 
frequency, multiple model frequencies are possible 
candidates (not necessarily only the closest one). 
This is particularly problematic in cases were no 
prior mode identification is available. 

3. As the observational uncertainties decrease, the 
contrast in x 2 (and even more so the contrast in 
probabilities) between different models increases. 
If the model that minimizes/maximizes Equa- 
tion ^ /Equation ^ is not the correct model due 
to missing physics, this increase in fitting contrast 
is misleading and unwarranted. 

4. In static grids the finite grid resolution increases 
the risk of missing the most adequate model that 
the code could produce. If there are systematic 
differences between even the most adequate mode 
and the observations, the "contrast enhancement 
effect" will be magnified. For the same reason, 
adaptive grids run into the same problem and will 
miss the correct parameter space region to finer re- 
solve in the first place. 

As a consequence of all these shortcomings, it is clear 
that a method is needed that considers the possibility of 
systematic differences. It is also mandatory to consider 
the finite resolution of our model grids. Solutions to these 
problems are presented in the following sections. 

3.2. The argument for probabilities 

There are obvious benefits to quantifying the best fit 
and the uncertainties in terms of probabilities. With 
probabilities for each specific model, we automatically 
obtain probability distributions for each parameter of the 
model itself. We can furthermore consistently compare 
different grids and see which set of input physics is more 
probable, given all our current information and the data. 

However, there are much stronger arguments for a 
probabilistic approach. Marginalization allows us to con- 
sistently treat nuisance parameters, while the sum and 
product rules allow us to clearly formulate the question 
we are asking. This question is "Given the observed 
frequencies, our knowledge about the star and model 
physics, which model(s) best represent the star in terms 
of its fundamental parameters and general physical prop- 
erties as probed by the pulsation modes?" In reality, this 
general question has to be further refined as we encounter 
more complicated situations like: "We have model fre- 
quencies that could potentially show negative or pos- 
itive systematic offsets, or no such offset at all, when 
compared to our observations. They could be influenced 



by rotation or actually be rotationally split frequencies 
themselves. They could be bumped I = 1 modes or I = 
modes. Given all of these possibilities, which model is 
the most adequate one, and how well is the solution con- 
strained?" . From the viewpoint of probability theory the 
only way to treat such a set of possibilities and get mean- 
ingful answers is to use the sum rule and product rule, 
as we will show in the next section. 

3.3. Ambiguous mode identification 

As a first improvement to the general approach of as- 
teroseismic model fitting, we can involve the sum rule to 
consistently consider uncertainties (or even ignorance) in 
our mode identification. In essence, if there is no unique 
proposition /i ;C »-H,m, Equation (|6| changes to 

lobs ( " mat c h ~j 

P(D\Mj,I)= H \ J2 P(ho,ho^k, m \Mj,I)\ (7) 
with 



i=l v. k=l 



P\fi,oi fi,oi-tk,m\Mj , I) — 
P(Ji ,oi— ^fe,m 

Mj,I). 



(8) 



Here the sum over the index k means that all possible 
and mutually exclusive assignments n matc h of one ob- 
served mode to a number of calculated frequencies fk m 
have to be taken into account as an "or" propositior 
Note that due to the product rule of probability, the 
terms in each sum now include the conditional prior prob- 
abilities P(fi, i->-k,m\Mj,I). These have to be normal- 
ized so that J2liT h P(fi,o^k,m\Mj,I) = 1. The most 
conservative assignment is to assign equal probabilities 
P(fi,oi->k,m\Mj,I) = l/n mat ch to each possible scenario. 
However, if more information is available (e.g., a mode 
could be identified to be either I = or I = 2 with specific 
probabilities for both cases as found by some peak bag- 
ging program), this can easily be encoded at this stage. 

The end result is a product of weighted sums of prob- 
abilities, where the weights are given by the respec- 
tive prior probabilitiesrj This product is the correctly 
normalized likelihood for obtaining the data, given the 
proposition that any one of the proposed scenarios is cor- 
rect. Note that if there is an unambiguous assignment 
/i,oi-n,m f° r every observed frequency, each prior proba- 
bility P(fi, h±Lm\Mj,I) — 1 and Equation afl simplifies 
to Equation ([6]) . Now that we have includea our uncer- 
tainties concerning the assignment of model frequencies 
and observed frequencies, we will deal with uncertainties 
in the validity of the model frequencies themselves in the 
next section. 

3.4. Treatment of systematic errors 

As a next step, we now show how to treat the prob- 
lem of imperfect models. As mentioned before, apply- 
ing standard techniques that rely on minimizing the 

3 Hereafter, a possibly ambiguous frequency assignment will al- 
ways be denoted as /i i->fc in- 

4 A common misconception is that these "priors" are only there 
to allow us to incorporate prior information. In reality, they are 
formally required by the product rule and ensure that the result of 
Equation |7| is always properly normalized. 
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quadratic differences between the observations and the 
models will give incorrect results if systematic differences 
exist. The alternative of correcting for such imperfec- 
tions prior to modeling is also undesirable if the correc- 
tion is not known to be universally applicable. 

To treat any systematic deviation from the model fre- 
quencies due to unmodeled physical effects, we simply 
expand the models Mj by considering an additional sys- 
tematic error parameter for each tested frequency. The 
aim is to construct new values Aa to compare with the 
observations according to 



fi,A = Am + l^i 



(9) 



Here, Aj is the absolute value of the systematic error. 
7 = for7 = — f and determines whether the model fre- 
quency is expected to be systematically higher or lower 
than the observed frequency. To keep our notation from 
occupying too much space, we will implicitly assume the 
value of 7 to be constant throughout the following deriva- 
tions, and attribute this to our prior information i\ Aj 
is an unknown parameter but as long as its lower and 
upper boundaries can be roughly estimated, it can be 
treated fully consistently in the probabilistic framework. 

In the following, we will again work out an example of 
only one observed and calculated frequency. Therefore, 
for the derivation the assignment Aoi->-i,m is unique. We 
will then provide the extension to multiple frequencies 
and ambiguous mode identifications. 

Using the new parameters, the equivalent to Equa- 
tion ^ is 



P(Ao,Aj|/j 

,01-4-2, m ) 



1 ,M 7 A ,/)P(/ i , |A i ,/, 



i, oi— H.in 



,Mf,I) 



P(A l |/ l , ^, m ,M. A ,/)x 



f 



2na. 



exp 



(/»,o ~ Am - l&iY 
2a? 



(10) 



Here the symbol M A simply denotes the model Mj 
augmented by the new parameter A, . Self-evidently, the 
product rule again requires that we introduce a prior 
probability P(A i \f io ^ im ,AIj i ,I). This can either en- 
code prior information about the expected behavior of 
the error, or be simply assigned by considerations of sym- 
metry. Again it is required that the integral over the 
prior / P(A i |Ao^,m,Mf ,J)rfAi = 1. 

It would now be possible to try to find the Aj that 



maximizes P(Ao, Ai\fi, ol -n, m , Mf, I) in Equation \l0\ 



However, this is completely irrelevant for our needs. 
In case of multiple observed frequencies it would also 
quickly lead to a highly dimensional parameter space 
that we are not interested in navigating. Instead, we are 
interested in finding the probabilities of the models M^. 
To do this it is necessary to integrate out Aj which we 
have just introduced. We obtain the marginal likelihood 

P{fi,o\fi,ot-^i,mjMj , I) = 

/•Ai.max (11) 

/ P(Ao,A i |/ i M?J)d&i. 



This integral naturally depends on the shape of the prior 
probability distribution for Aj, and can easily be evalu- 
ated numericalljr] It represents the likelihood of obtain- 
ing the value of the observed frequency Ao given that Mj 
predicts a frequency Am but that there is a possibility of 
a systematic difference Aj, between Aj jm ; n and Aj max . 
Furthermore, it is fundamentally constrained and prop- 
erly weighted by the prior we assigned. This result is 
now easily extended to multiple modes and ambiguous 
mode identifications. Equation Q becomes 



™obs f "match 



P(D\M*,I) = 1] 1 E P (Ao, Ao^fe,m|Mf ,1) 



i=l (. k=l 



and 



P(/i,oi fi,o^-k,m\Mj , I) — 

P(Ao^fc,m|M/,7)P(Ao| Ao^,m, Mf , I). 



(12) 



(13) 



In summary, we have to calculate a prod uct of weighted 
sums of integrals in the form of Equation (111, where the 
summation is performed over every possible assignment 

fi,ot- >fe,m* 

Note that even when we choose to consider systematic 
deviations, we usually do not expect them to be signifi- 
cant for all frequencies. For good models some frequen- 
cies should already match well "right out of the box" . In 
particular, this is true for all frequencies in the idealized 
case where we have (finally) found a way to correctly 
model all the effects that previously caused systematic 
deviations. 

One might think that this is taken care of by 
setting Aj an i n = 0. However, unless the prior 
P(Aj| Ao>->.fc,m, Mj, I) is a <5 function at Aj = 0, it is 
much more likely that Aj > 0. This means that a priori 
a model will be preferred which shows at least a small 
deviation from the observations, depending on the ob- 
servational uncertainties and the steepness of the prior. 
The limiting case however, the 8 function, corresponds 
to a whole different model which is simply the standard 
model without systematic deviations, Mj. Thanks to 
the sum rule, there is an elegant solution for taking this 
alternative into account. 

For the mutually exclusive logical propositions^] M A 
and Mj we can calculate 

P(Ao: fi,o^k,m\M^ + Mj,I) = 

\I) + P{M j J it0 ,f i 

,oi— m 



P{Mf\I) + P{Mj\I) 



P{M?\I) 



P(Mf\I) + P(Mj\I) 



P(Ao,Ao^fe,m|M„ A ,/) + 



P(Mj\I) 



P(M A \I) + P(M 3 \I) P{k °' ko^k,m\Mj,I). 



(14) 



5 For several simple shapes, such as the beta prior introduced in 
the next section, there also exist analytical solutions. 

6 Note that from a logical standpoint even if A = 0, Mf- is still 
a different model than Mj because the prior is not a <5 function. 
Therefore, they are always mutually exclusive. 
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Note that here Mf + Mj means "M^ or M i is true" . 
This is the likelihood of observing the frequency value 
fi o) given that a systematic deviation either does or does 
not exist. The principle of indifference as the most con- 
servative approach for the prior probabilities obviously 
demands P(Mf-\I) = P(Mj\I) = 0.5, but if more infor- 
mation is available, it can be encoded here. This result is 
also easily generalized to the case of multiple frequencies 
and ambiguous mode identification. 

3.5. The choice of the prior for Aj 

A very important detail to consider when extending the 
models with systematic error parameters is their prior 
probabilities P(Aj|/j i0h ^.ji iro , 7). There is a basic 
choice between two possibilities. The first is to use un- 
informative (or ignorance) priors, or alternatively, maxi- 
mum entropy priors. Uninformative priors can be derived 
from arguments of invariance to specific transformations, 
while maximum entropy priors should satisfy the maxi- 
mum entropy criterion for a given set of constraints. The 
other possibility is to use priors derived from heuristic or 
physical arguments. 

The specific form of the prior probabilities of Aj are 
part of the model that is evaluated, as indicated by 
the notation^ They are not necessarily "prior" as in a 
sense of "before obtaining observations" , but conditional 
probabilities required for the correct normalization, as 
demanded by the product rule of probabilities. They 
encode specific ways in which we expect Aj to behave, 
given our grid of frequencies and our information (which 
of course can be influenced by previous observations). 
For instance, if we expect our best model to minimize 
the systematic deviations, the prior should assign larger 
probability densities to smaller Aj, so that models with 
smaller deviations will be more probable. On the other 
hand, if we expect our best model to show more erratic 
deviations, a flat uninformative prior is a better choice. 
After a complete evaluation of the probabilities and like- 
lihoods, the Bayesian evidence will indicate whether the 
state of information encoded by the priors is supported 
by the data or not. 

As a first important example of an uninformative prior, 
consider a uniform prior 



P(&i\fi, 



OH >k,m i 



Mf,7) 



A ? ; 



= const. 



(15) 

This means that all values of Aj are equally likely. With 
such a prior, every model that predicts frequencies at any 
value between fi l0 + Aj jm i n and /j, +Ai )max has the same 
maximum lik elihood (i.e., the same maximum value for 
Equation ([To])). 

On the other hand, a Jeffreys prior 



P(A i |/ i , ^, m ,Mf,7) = 



In (Aj iiriax / Ai,min) 



, (16) 



assigns equal probability per decade and, in terms of the 
probability density, favors smaller values of Aj. This 
prior is obviously not defined for Aj — 0, i.e., it requires 



7 The prior is described by P(A;|/; OH 
e.g., P(Ai\I) 



k,mi /) rather than, 



Aj^min > 0. This is problematic for, e.g., surface ef- 
fects that approach zero at low orders. However, when 

Aj. m i n — 

given by 



Aj im i n = 0, one can use a modified version of this prior 



P(X\fi. 



Mf,I) = 



1 



(A, + c)ln[(A 4 



+ c)/c]' 

where c is a small constant. For values smaller than c, 
this prior acts more or less like a uniform prior, while 
for higher values it behaves like the usual Jeffreys prior. 
This prior is n owadays often used in "peak-bagging" al- 
gorithms (e.g., [Gruberbauer et al. 2009; Benomar et al. 



2009 



no 



»• 

ob 



Handberg & Campantc 20rnH However, there is 
jective criterion for how to set c, and various tests 
we conducted with our grid fitting code have shown that 
the choice of c can have a large impact on the evidence 
values. 

Consequently, we argue that any priors used for a sys- 
tematic error parameter Aj = [0, Aj. max ] should be func- 
tions that are clearly defined by the parameter limits, 
without additional parameters that have large effects on 
the evidence. The uniform prioi[^] is such a prior, as are 
priors derived from the beta distribution (given in units 
of our problem) 



7 J (Aj|/ l 

,oi— >/c,m: 

M A , 7) cx 



A, 



(18) 

With a = 1 and j3 = 2 this simply leads to a linearly 
decreasing probability density 



P(Ai\fi t0h + kim ,Mf,I) 



V2 



(A, 



A,. 



(19) 

It is the only prior that allows for a linearly decreasing 
probability density, is properly normalized, and reaches 
zero at Aj = Aj jmax . It also lead s to an analytical solu- 
tion for the integral in Equation ( 11 1. Thus, it satisfies 
all our requirements for a prior with which to minimize 
systematic errors. 

W e ha ve compared the results obtained from Equa- 
tion ( |19[ ) (hereafter: beta prior) with several other plau- 
sible choices, such as an exponential distribution with 
expectation value Aj )max /2 and a modified Jeffreys prior 
with c = CTj m , and we find them to yield comparable 
results and evidence values. Due to the clarity of its def- 
inition and lack of additional parameters, we therefore 
argue that the beta prior is an appropriate choice for a 
non-flat prior. We will show how to use it in Section[2]in 
a worked example. 

Lastly, priors based on heuristic or physical arguments 
obviously vary strongly with the specific problem to 
which the fitting method is applied. As an example, 
when modeling surface effects on p-mode frequencies, the 
prior could be Gaussian , following the heuristic frequency 
correction proposed by Kjeldsen et al. (2008). It would 
be a function of frequency, expecting greater deviations 
toward higher-order modes. The width of the Gaussian, 

8 In fact, the uniform prior is consistent with a beta distribution 
with a = 8 = 1. 
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however, would be again a rather arbitrary choice, lead- 
ing to potentially different evidence values. Such priors 
clearly need to have a strong basis either in theory or 
prior observations. 

3.6. Bumped modes and finite grid resolution 



Equation ( 12 ) represents the final likelihood for obtain- 
ing the observed frequencies given our (extended) model 
Mj . This model still represents only a single point 
in a discrete gridr] However, the probability is small 
that a single model in the grid corresponds to the "true 
best model" our code can produce. The problem be- 
comes worse as the grid resolution is lowered, or as mode 
frequencies are changing quickly or unpredictably from 
one model to the next (e.g., avoided crossings, magnetic 
shifts). The probabilities (or x 2 -vahies) we obtain will 
not be a fair assessment of the model physics, even at 
higher grid resolutions. Even worse, the overall evidence 
for the grid will be finely tuned to the positions of all 
models in the grid. This makes it difficult to compare 
different grids with different physics. We will now show 
how to improve on this. 

In a sequence of models along a single evolutionary 
track, except for the first and last models, each model 
Mj has two neighboring models Mj-\ and Mj+i. In 
most cases these adjacent models will contain the same 
modes, and their changing values can be traced from 
Mj_i to Mj and Mj+i. Now we declare the difference 
between observed and calculated frequency as a new free 
parameter 

Sfi = fa - fa- (20) 

This value is fixed if only a single grid point is consid- 
ered. However, we can split the evolutionary tracks into 
segments in between grid points, and define 



3fi,j— fi,o 

and equivalently 

$fa+ = fa 



fi,Mj-i + fi,Mj 



(21) 



(22) 



Adding Sfi as a new parameter to the equations de- 
rived in the earlier sections, we change our focus to 
evaluate probabilities of model track segments T^ A cen- 
tered around the models M- . To do this, we again use 
marginalization to integrate out both Aj and Sfi to ob- 
tain the marginal likelihood. We obtain 



P(fi,o\fi,o\ 

fA 4 . max r 6f i:j+ 

Sfi. 



T) 



P(fa, A», <y/i|/i,oi->-t,m, T } , I) dAi dSfi. 

(23) 



If the priors for Aj do not vary greatly from one model 
to the next, Aj and Sfi can be considered to be inde- 
pendent parameters. It is therefore possible to use the 
product rule to separate the conditional probabilities 

9 Note that this is also the case for approaches using an adap- 
tive grid, since each iteration of an adaptive scheme is based on a 
discrete grid. 



^(AojAj, | /i,oi->i,m j Tj , /) 
P{^-i\fi,oi-^i,mj Tj ,I)x 

P{Sfi\fi^i, m ,Tf,I)x 

-f (/i,o| Aj, Sfi, /i,oi->i,m! T) i I) 



(24) 



Furthermore, since we evaluate the complete evolu- 
tionary track segment we can assume a uniform prior 
probability P(Sfi\fa^ itm ,Tf,I) = 1/ (Sfa + - 8 fa-). 
With these definitions, the integral over Sfi can easily be 
calculated analytically. The equivalent to Equation ( |To| 
becomes 



P{fi,oiAi\fi oi -^i m: Tj , /) 

P{^i\fi,oi-ti,mi Tj , I)X 
1 

V 

2 i S fa+ - 6 fiJ-) 



(25) 



erf 



&fa+ - 7^ 



rf 



/ 8 fa- - 7A, 



where erf is the error function. The remaining integral 
over Aj again has to be carried out numerically. Figure|2] 
shows an example for the definitions introduced above, 
given three models in a solar evolutionary track. 
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Figure 2. Example for the definition of 5i j— and 5j j-f. (see the 
text). Four radial orders of I = 1 modes from three adjacent models 
in a high-resolution grid of solar models are shown. The triangles 
represent the central model. The frequencies for the adjacent mod- 
els along the evolutionary track sequence are depicted as squares 
and white circles. Black circles (and error bars) indicate observed 
frequencies published in |Broomhall et al.| ( [2009^ . 8ij— and <5i,j+ 
for a single mode are represented using arrows. The insert shows 
an unzoomed version of the I = 1 and I = 3 ridge. 
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We have now used the free parameter Sfi to "trace" 
each mode through segments of the evolutionary track, 
and compare it to the observed frequencies, retaining 
the possibility of systematic differences. Note that our 
only assumption here is that the mode frequencies change 
smoothly between the frequencies given by the constrain- 
ing models. In principle, this approach can be carried 
out in multiple dimensions (e.g., not only along the evo- 
lutionary track in stellar age but also between tracks in 
mass). As before, an extension to multiple frequencies 
and ambiguous mode identifications is straightforward. 

We stress that this approach only locates the region of 
highest probability given the current grid, and given un- 
specified behavior of frequencies in between grid points. 
It is thus best used for frequencies whose behavior is dif- 
ficult to capture, e.g., due to mode bumping, or for a 
first general assessment of a very coarse grid. Given a 
dense enough grid, regular frequencies that are expected 
to change approximately linearly from one grid point to 
the next need to be treated using interpolation, since the 
integration over the model gaps for individual modes, in- 
dependently of all other modes, would allow for highly 
unphysical models. 

Therefore, in order to obtain a final best model and 
uncertainties for the model parameters, the regions of 
substantial probability should be further refined after 
the track probabilities have been calculated. Eventu- 
ally, the grid is resolved enough so that well-defined dis- 
tributions arise. In dense enough grids, this can easily 
be accomplished by interpolation of the frequencies in 
between grid points without violating the condition of 
hydrostatic equilibrium. This can also be done during 
run-time with arbitrary precision using probabilities, by 
interpolating between grid points and using the sum rule 
to calculate a probability representative of the original 
grid resolution using the interpolated models. Naturally, 
modes that change erratically, should be excluded from 
such an interpolation routine and treated as shown above 
instead. 

3.7. Model probabilities 

So far we have only shown how to calculate the likeli- 
hood for standard pulsation models, models that contain 
systematic differences, and also evolutionary track seg- 
ments. In order to obtain the probabilities for individ- 
ual models (or track segments), we want to use Bayes' 
theorem, assign model priors, and calculate the total evi- 
dence for each model grid. The simplest method to assign 
model priors in the absence of any other prior informa- 
tion, is to use the principle of equipartition and assign a 
uniform prior 



density for the model temperature is 



P(M,-|i) = 1.0/iV M , 



(26) 



where Nm is the number of models (or, equivalently, iVr 
would be the number of track segments) that are ana- 
lyzed. 

Although each model or track only predicts a number 
of frequencies, it implicitly represents values or ranges for 
fundamental parameters like T e g or L, which can be com- 
pared to (and constrained by) different and non-seismic 
observations. For instance, assuming our prior photo- 
metric and spectroscopic observations of a pulsating star 
indicate T c s = T spcc ± a spec then the prior probability 



P(T oB , j \I) = kexp 



T eS,jY 



2<ec 



(27) 



This example assumes that the uncertainty in T c g follows 
a Gaussian distribution. A: is a normalization constant 
depending on the absolute lower and upper plausible lim- 
its Of T e ff. 

All the different implicit parameters for which prior 
observations or other fundamental constraints are avail- 
able, and hence prior knowledge exists, then can be used 
for prior probabilities which combine into an overall prior 
for model Mj. As an example 



P{M 3 \I) = P(T effj |/)P(L#)P([Fe/Hl . \I). 



(28) 



if we assume separable priors for simplicity. If proba- 
bilities of track segments are calculated, such a prior 
could be approximated by a product of separable inte- 
grals, which are easily evaluated analytically. Again, in 
order to obtain a proper prior and therefore proper val- 
ues for the evidence, the integral of the prior probability 
over all possible models/tracks in a grid should be 1. 
By calculating, e.g., 



P(Mj\D, I) 



P{Mf\I)P{D\Mf,I) 



(29) 



or, e.g., in the case of rapidly changing modes with or 
without systematic errors, 



P(Tf+Tj\D,I) 



P(TA + T 3 \I)P(D\T : A + Tj,I) 

Ek=i P (T£ + T 3 \I)P(D\T* + T~I) 

(30) 



we obtain the probability of M- or, respectively, T^+Tj 
given our prior knowledge (or lack thereof) , our grid, and 
the set of observed frequencies. Note that the denomi- 
nators of these equations represent the evidence or likeli- 
hood for the grid as a whole. We can therefore use these 
as likelihoods when we want to compare different grids 
with different input physics. 

4. APPLICATION TO SURFACE EFFECTS 

As mentioned in the introduction, shortcomings in 
modeling the outer stellar layers produce systematic de- 
viations in comparison to the observations. These devi- 
ations seem to be such that model frequencies tend to 
be higher than the observe d frequencies, and the refore 
7 = -1 (see Equation {§])). |Kjeldsen etaL] ( |2008[) have 



proposed to calibrate a power-law description of the de 
viations by measuring the surface effects in the Sun, and 
then fitting this relation to frequencies of other stars. 
Their correction expressed in terms of our definitions has 
the general form of 



7A* 



fi , m 
/ref 



(31) 



where / re f is some reference frequency, typically the fre- 
quency of maximum power i/ max , and a and b are pa- 
rameters to be fitted. From their fits, |Kjeldsen et al.| 
determined b w 4.90 in the Sun, which has subsequently 
been used for other stars by a number of authors. A 
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comprehensive implementation of this formalism into a 
Y 2 -fitting algor i thm w as presented in a recent study by 
|Brandao et al.| ( |2011[ ). However, even in this more ad- 
vanced approach, there is still a choice of a and b re- 
quired. Moreover, complications for modes of different 
spherical degree and also bumped modes arise because 
they do not necessarily conform to this relation. The 
authors propose to alleviate these problems by introduc- 
ing additional model-dependent parameters that approx- 
imately correct for some of these deviations. While this 
approach is a great improvement over applying a fixed 
surface-effect correction (or no correction at all), our ap- 
proach is much more powerful. It allows for a much 
greater flexibility and leads to clearly defined probabilis- 
tic results. 

4.1. Priors for surface effects 

As we want to treat systematic errors of more or less 
unknown magnitude, the most ge nera l approach is to 
use the flat uniform prior (Equation ( 15 1). Imposing onh 



minor additional constrai nts, as we argued in Section |3.5[ 
the beta prior (Equation ( |19[ )) can also be used to give 
more weight to models which minimize these unknown 
errors. We can use both priors and compare the Bayesian 
evidence to tell us which interpretation of the surface 
effect is better supported by the data, given our model 
and everything we know. Moreover, irrespective of which 
prior is chosen, we also always allow for the pos sibil ity 
of no surface effects at all, as discussed in Section [3~4| 

This now gives us enough flexibility to consider a pos- 
sibly frequency-dependent behavior of the surface effects. 
However, instead of "predicting" the behavior of Aj as 
is done by modeling the surface deviations through a 
power law, we will prescribe the behavior of its upper 
limit A ijmax . In contrast, the lower limit should always 
remain 0, since our ultimate goal is to find models that 
correctly describe the surface layers (and therefore ap- 
proach Aj — » 0). 

The choice of the largest allowed A max is not unique, 
but it should be sensible and used consistently through- 
out the analysis. A reasonable strategy is to use 
sup(A max ) = Av, the large frequency separation of each 
specific model, as a sensible upper limit. If the systematic 
differences between observations and models are larger 
than the average distance between modes of adjacent ra- 
dial order, we no longer recognize this as a valid fre- 
quency assignment^! With this upper limit defined, we 
now want to modeTdifferent types of surface effects. If 
we have no preference for any frequency-dependent trend 
(i.e., all we know is that observed frequencies are lower 
than model frequencies) we require that all frequencies 
have equal A max . 

On the other hand we ca n also use a more specific 
model, such as Equation (31 1, but retain the sa me fl exi- 



bility. The surface effect as shown in Equation (31) de- 
pends on two parameters. The power-law exponent b 
determines how quickly the surface effect increases as 
we move to higher frequencies, whereas a is simply a 
scaling factor. We are not interested in the scaling pa- 
rameter, since the scaling (i.e., the magnitude of the off- 
set) is governed by our condition that for each model 
sup(A max ) = Av. It is taken care of by the fact that we 

10 This condition may be relaxed at the highest radial orders. 



are marginalizing over A^ anyway. Since the largest sur- 
face effects are expected at the highest frequency / maXjm 
in the model, it follows that for a specific b 



fi.i 



(32) 



Figure[3] shows how these definitions affect the prior prob- 
ability density P(A i |/i i0 ^i jm , Mj A , /) as we increase the 
value of Ai for both the constant and the power-law ap- 
proach. With all the A maXj j set, we can then use the pri- 
ors as discussed above for all our calculations. Note that 
we can also easily evaluate new composite propositions 
at this stage and compute the probability for a hypoth- 
esis that allows for, e.g., a range b = {4.4,4.65,4.9, ...}. 
This is done i n th e same way as was explained earlier 
(see Equation ( 14 1). 



4.2. Detailed analysis of the Sun 

As an example for how to implement the surface-effect 
treatment, we will consider the solar I = 0, 1, 2, and 3 p 
modes obtained by using BiSON data ( |Broomhall et al.| 
2009 1 . For our models, we used a large and dense solar 
gfidobtained with YREC (|Demarque et al.| [2008). The 
model grid spans: masses from 0.95 Mq to 1.05 Mq in 
steps of 0.005 Mq, initial hydrogen mass fractions from 
0.68 to 0.74 in steps of 0.01, initial metal mass fractions 
from 0.016 to 0.022 in steps of 0.001, and mixing length 
parameters from 1.8 to 2.4 in steps of 0.1. These param- 
eters are also summarized in Tablc[T] 

Our model trac ks begin as completely convective L ane- 
Emden spheres ( |Lane||l869[ jChandr asckhar 195 7]) and 



are evolved from the Hayashi track flHayashil 1961 1 



through the zero-age main sequence (ZAMS) to 6 Gyr 
with each track consisting of approximately 2500 mod- 
els. Only models between 4.0 and 6.0 Gyr are in- 
cluded in the model grid. Constituti ve p hysics include 
the OPAL 98 j|glesias fc Rogers||1996[ ) and | Alexander fc 
Ferguson] dl994[) opacity t ables using the GS98 mixture 
i Grevesse &: Sauval||1998[ ), an d the Lawrence Livermore 
2005 equation ot state tables ( Rogers|1986 Rogers et al. 
19961. Convective energy transport was modeled using 
the B ohm-Vitense mixing-length theory (Bohm-Vitense 
19581. The atmosphere mode l follows the ('J'-t) rela- 
tion by Krishna Swamy l (|1966 ). Nuclea r reaction cross- 
sections are~Eom~^^^^^^^. 2001). The effects of 
helium and heavy element diffusion ( Ba 
were included. Note that our atmosphere models and dit 
fusion effects have been shown to require a larger value of 

mixing length parameter (a m i ss 2.0 2.2) tha n stan- 

dard Eddi ngton atmospheres (a m i « 1.7 1.8) ( |Guen 

ther et al ||1993[ ) 



1'he pulsation spe ctra were comput ed using the stel- 
lar pulsation code of Guenther (1994), which solves the 
linearized, non-radial, non-adiabatic pulsation equations 
using the Henyey relaxation method. The non-adiabatic 
solutions include radiative energy gains and losses but 
do not include the effects of convection. We estimate the 
random la uncertainties of our model frequencies to be 
of the order of 0.1 /LtHz. 

We analyzed our grid using adiabatic and non- 
adiabatic frequencies, and employed three different 
surface-effect models: 
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Figure 3. Behavior of the beta prior for systematic offsets in an cchcllc diagram. The squares represent model frequencies, while the 
shaded "trails" indicate the prior probability densit y for varying A;. The left panel uses equal A max , whereas the right panel uses a 
power-law A max with exponent b = 4.9 (see Section [4.1[ |. Note that the uniform prior is not shown, since it simply assigns a constant 
probability density. 



Table 1 

Parameter Ranges for the Solar Grid. 



Parameter 


Range 


Step Size 


Mass 


0.95-1.05 


0.005 


X 


0.68-0.74 


0.01 


Z 


0.016-0.022 


0.001 


"ml 


1.8-2.4 


0.1 



(c) same as (b) but with an additional Gaussian 
constraint on the age of 4.603 ±0.0075 Gyr, derived 
fro m the estimated age of th e solar system found 



by Bouvier & Wadhwa 
main sequence phase o: 



(2010) and an average pre- 



Notes. Masses are given in units of solar masses; a m i is the 
mixing length parameter. 



Ml: frequency- independent surface effects 



with 



M2: frequency-dependent, "canonical" surface ef- 



fects with Aj max following Equation ( 32 ) with b = 
4.90 

• M3: same as M2, but with b as a free parameter 
marginalized from b = 3.0 to b = 6.0 

For each frequency evaluated throughout our model 
grid, irrespective of the surface-effect model, we also con- 
sidered the possibility of no surface effect, i.e., we consis- 
tently calculated P(Mj+Mj\D, I). To take into account 
the discrete nature of the grid, we interpolated along the 
evolutionary tracks during run-time by a factor of 20, 
thereby increasing the effective "frequency resolution" of 
the grid to below the random uncertainties of the model 
frequencies. All models were evaluated with 

• (a) a uniform prior for all track segments 

• (b) a prior using normal distributions for the 
constraints M = 1.0000 ± 0.0002 M , logT cff = 
3.7617±0.01, andlog(i/i ) = 0.00±0.01, where 
YREC uses the followi ng adopted values for M n = 
1.9891 ± 0.0004 • 10 33 g ([Cohen fc Taylor|l986| and 
L = S.SSlSiO.OOg-lO^ergs" 1 (the average of the 
ERB-Nimbus and SMM /ARCRIM measurements; 
Hickey fc Alton] ( [l983| ) 



our models of 35 ± 5 Myr. 

For the A, we consistently used beta priors, as dis- 
cussed in the previous section. Our calculations yield 
the most probable models and uncertainties for all these 
approaches, and they also give the Bayesian evidence for 
each approach. The results are summarized in Table[2j 
We also computed the probabilities using uniform priors, 
but found similar results with lower evidence (several or- 
ders of magnitude) values than for the corresponding beta 
prior analysis. 



Table 2 

Evidence for the Solar Grid Using 



the BiSON Data Set 



Surface 


HRD 


log 10 {evidence) 


log 10 {evidence) 


Model 


Prior 


(adiabatic) 


(non-adiabatic) 


Ml 


a 


-233.4 


-229.9 


M2 


a 


-189.8 


-186.7 


M3 


a 


-189.8 


-187.2 


Ml 


b 


-235.5 


-231.6 


M2 


b 


-190.8 


-187.1 


M3 


b 


-191.4 


-187.9 


Ml 


c 


-236.7 


-235.1 


M2 


c 


-193.5 


-190.7 


M3 


c 


-192.6 


-189.4 



Notes. Sec the text for the definition of models and prior a, 
b, and c. Results for models M2a, M2b, and M2c, which are 
analyzed in more detail, are indicated in bold face. Note that 
small numbers are expected. 

The non-adiabatic frequencies consistently produce 
larger evidence values than for the respective adiabatic 
case. This is no surprise, as the non-adiabatic frequen- 
cies are in general better at reproducing the higher fre- 
quencies. Overall, model M2a shows the largest evi- 
dence, followed by M2b and M3a. Note that Mia, 



Toward a New Kind of Asteroseismic Grid Fitting 



11 



Table 3 

Most Probable Models for the Complete BiSON Data Set Model Fitting. 



Model 


Mass 


Age 


X 


Zq 


z s 


"ml 


Probability 


M2a 


1.015 


4.885 ± 0.006 


0.73 


0.017 


0.0153 


2.2 


0.54 




1.005 


4.713 ±0.006 


0.72 


0.017 


0.0153 


2.2 


0.21 


M2a 


1.000 


5.017 ± 0.006 


0.72 


0.018 


0.0161 


2.1 


0.68 




1.000 


4.983 ± 0.006 


0.71 


0.019 


0.0160 


2.2 


0.17 


M2a 


1.000 


4.591 ± 0.005 


0.72 


0.016 


0.0144 


2.2 


0.95 




1.000 


4.562 ±0.005 


0.71 


0.017 


0.0153 


2.3 


0.05 



Notes. 

years less 

with respect to the specific surface-effect model and prior combination 



Age is given in billion years and is computed from the pre-main-sequence birthline. The age from the ZAMS is 35 ± 5 million 
Xq and Zq are initial hydrogen and metal mass fractions, Z B is the metal mass fraction in the envelope. Probabilities are given 



Mlb, and Mlc, which use frequency-independent pri- 
ors for the surface effects, and therefore are extremely 
flexible, fail compared to the other models. Also, M3a 
and M3b cannot beat their M2 counterparts. These are 
examples of how marginalization and the consistent nor- 
malization of probabilities work together to penalize more 
flexible models if they cannot generate considerably bet- 
ter results. Model M3c has a greater evidence than M2c, 
but the most probable stellar models are the same in 
both cases, suggesting that these models fit well, but do 
not necessarily adhere in detail to the standard surface 
correction. 

At first glance it might be unsettling that M2a has 
a slightly greater evidence than M2b (and significantly 
greater evidence than M2c). This indicates that there 
are models in our grid which reproduce the pulsation 
spectrum very well but do not match the solar funda- 
mental parameters. A correctly calibrated grid would 
produce higher evidences with a prior restricted to the 
true solution. However, regardless of whether or not 
we include the fundamental parameter constraints, we 
are still finding models that match the oscillations con- 
straints reasonably well. Furthermore, recall that the ev- 
idence is only the likelihood of obtaining the data, given 
that the approach is correct j^j] We know that the a prior 
is misrepresenting our state of information. The solar 
prior approach b more correctly encodes what we know 
about the Sun, and the age prior c puts even tighter 
constraints on the pulsation models. Ignoring this in- 
formation (using prior a and setting equal conditional 
probabilities) is an interesting and necessary exercise to 
study the consistency of the results, and how the differ- 
ent models, approaches, and priors work. For an actual 
detailed study of the solar model physics, however, it is 
not appropriate. We can nonetheless compare the re- 
sults, restricting ourselves to the non-adiabatic frequen- 
cies and the on average best model for each prior, M2. 
The resulting parameters are displayed in Table[3] 

The results obtained without using our prior Knowl- 
edge of the Sun for model M2a are spread over several 
models in the parameter space that can fit the observa- 
tions quite well. However, for the most probable models, 
the mass and age are inconsistent with our prior knowl- 
edge. These models seem to produce smaller surface ef- 
fects, and are therefore preferred. For model M2b the 

11 In order to obtain correctly normalized probabilities for the 
different approaches themselves, we have to introduce conditional 
probabilities like, e.g., P(M2a|7) or P(M2b|I) and use Bayes' 
theorem. Only comparing the evidence amounts to setting these 
conditional probabilities to be equal for all tested hypotheses (e.g., 
P(M2a|J) = P(M2b|/)). 



situation is similar. Although the mass is now fixed to 
the true solar value, we do not obtain models that are 
consistent with the solar age. 

For M3c a single combination of physical parameters 
dominates the results and manages to fit well all the con- 
straints we impose (mass, luminosity, T e fj, pulsation fre- 
quencies, and age). Loosening the conditions on T e g and 
the luminosity does not significantly change the result. 
We have also tested slight variations of up to 20 million 
years in the age prior and do not find the result to be 
affected. In all cases, we recover a tightly constrained 
most probable model with Z — 0.016 and Z s = 0.0144, 
and an age of 4.591 ± 0.005 Gyr. We th erefore find a re- 
sult similar to Houdek & Gough (2011). Given that our 
models take 35 ± 51Vlyr to reach the main sequence, our 
result is also consistent with meteoritic age determina- 
tions of t he solar system to within s everal million years 
(see, e.g. Bouvier & Wadhwa 20101. However, we also 
recover Xq = 0.72, which leads to an initial helium mass 
fraction of Yq = 0.264(1). This is different c ompared to 
the valu e of Y — 0.250(1) that was found by Ho udek fc 

(|2009|). 

the 



Gough but more consistent with Asplund et al 
fitting 



observations to the adiabatic frequencies, 
including the age prior, we also recover the exact same 
model. We also tested how sensitive the grid is to the 
prior constraints in order to estimate the actual impact of 
the pulsation frequencies on the probabilities. If we only 
evaluate the combined priors, ignoring the frequencies 
but including the prior on the age, we obtain Xq — 0.71± 
0.01, Z = 0.019 ± 0.002, Z s = 0.017 ± 0.002, age = 
4.603 ± 0.008 Gyr, and a ml = 2.1 ± 0.2. This leads us 
to conclude that the frequencies have a decisive impact 
and actually select the low-metallicity models, no matter 
whether adiabatic or non-adiabatic model frequencies are 
used. 

However, it has to be stressed again that the evidence 
drops by almost two orders of magnitude when we intro- 
duce the age prior. This can be understood by the fact 
that the solution is so well constrained and at the edge 
of our current parameter space in Zq, and that many 
other models can also produce similar pulsation spectra. 
It could also suggest that we might not have covered 
the true best model parameters yet in our current grid. 
Therefore, our next goal will be to extend the grid to 
lower metallicities, and also include different abundance 
mixtures, but this is beyond the scope of this paper. 

FigureH compares the BiSON observations with our 
most probable model at the correct solar age. Even 
with non-adiabatic frequencies, significant surface effects 
can still be found. The measured surface effects them- 
selves are shown in Figure[5] together with least-squares 
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Figure 4. Non-adiabatic (shaded symbols) and adiabatic (open 
symbols) frequencies of the most probable solar model from eval- 
uating the BiSON frequencies (black circles + error bars) using 
approach M2c. Note that frequencies have been shifted upward 
by 5 fiHz, before calculating the x-axis values in order to prevent 
the 1 = 2 modes from wrapping around. 



fits fo llowing the relation proposed by |Kjeldsen et al. 
(2008). The magnitude of the surface deviations depends 
on whether the non-adiabatic or the adiabatic frequencies 
are used for the fit. Nonetheless, our method manages to 
identify the same exact model to be the most probable, 
even using the same surface-effect model, thanks to the 
power of marginalization. However, the non-adiabatic 
models are vastly preferred in terms of the Bayesian ev- 
idence. This is an example for how the presented ap- 
proach can be used to iterate toward improved stellar 
model physics, while still recovering meaningful stellar 
parameters from current asteroseismic investigations. 

We also determined surface-correction power-law ex- 
ponents for every spherical degree via least-squares fits. 
For both the non-adiabatic and adiabatic frequencies the 
best fitting exponents are mar kedly different from b = 4.9 
which was both advocated by |Kjeldsen et al. ( |2008[ ) and 
also used as the basis for our probabilistic surface model 
M2. This is also the reason why the M3c models have a 
greater evidence than their M2c counterparts. The fit- 
ted values range from b — 4.23 for non-adiabatic (I = 0) 
frequencies to b = 5.13 for adiabatic (I — 3) frequencies. 
Moreover, the power-law fits do not match the deviations 
very well at intermediate radial orders near 2400 fiHz. 
From our point of view, fixing the exponent t o b = 4.9 
for a least-s quares fit, as for instance done by Brandao 
et al. ( 2011 ), is therefore a potential problem since it does 



not even match the Sun very well, in particular when 
improved (e.g., non-adiabatic) physics are implemented. 
The probabilistic procedure has no problem with these 
deviations, even though it formally assumes an exponent 
of b = 4.9, since the magnitude of the surface effects is 
marginalized for every frequency. 
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Figure 5. Measured surface effects for non-adiabatic (filled cir- 
cles) and adiabatic (open circles) frequencies of the most proba- 
ble solar model from evaluating the BiSON frequencies using ap- 
proach M2c. The uncertainties of the differences are smal ler t han 
the symbols. Least-squares power-law fits (see Equation l |31| ) to 
the surface effects for the adiabatic (solid line) and non-adiabatic 
(dashed line) frequencies are also shown. 

4.3. Asteroseismic analysis of a Sun-like star 

To investigate the applicability of our method to cur- 
rent asteroseismic investigations, we also performed an 
"asteroseismic" analysis of a Sun-like star, simulated by 
artificially "degrading" the set of observed BiSON fre- 
quencies to a precision and accuracy expected from cur- 
rent space-based missions for average Sun-like solar-type 
oscillators. We first multiplied the uncertainties of the 
BiSON observations by a factor of 20, and then added 
corresponding random errors to the frequency values. 
Furthermore, we did not assume to have detailed prior 
information on the fundamental parameters. Instead, 
we fitted the "degraded" data set with a completely flat 
prior to the same grid as before, again using our surface 
effect model M2. 

Although a different most probable model is iden- 
tified, the overall results are comparable to our find- 
ings for M2a. They show a slightly larger spread of 
the model probabilities across the grid. Summarizing 
the uncertainties for the main parameters by calculat- 
ing the first and second central moments of the proba- 
bility distribution in our grid we approximately obtain 
M = 1.015 ± 0.007 Mq, age = 4.76 ± 0.10 Gyr, X = 
0.72 ± 0.01, Z = 0.017 ± 0.001, Z s = 0.0148 ± 0.0005, 
and a m i = 2.3 ± 0.1. 

However, these results become worse if we systemati- 
cally remove lower order modes which are crucial to "an- 
chor" the surface effect relation. To illustrate, we further 
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degraded our data set by only keeping 13 I = modes 
from 1950 to 3580 ^Hz, 12 I = 1 modes from 2020 to 
3505 /zHz, 10 I = 2 modes between 2080 and 3300 /iHz, 
and 8 Z = 3 modes from 2270 to 3220 /iHz. Similar data 
sets from Kepler and CoRoT with comparable uncertain- 
ties and numbers of modes have recently been analyzed 
in the literature. The results for the model parameters 
become M = 1.046 ± 0.007 M , age = 4.80 ± 0.43 Gyr, 
X = 0.72 ± 0.01, Z = 0.021 ± 0.01, Z s = 0.019 ± 0.001, 
and a m \ = 2.3 ± 0.1. Although the values are still within 
~ 5% we are almost at the border of our parameter 
space, and higher-mass models systematically outper- 
form lower-mass models. 

We know from investigating the BiSON data using our 
grid that we require a m i = 2.2 to fit all solar observ- 
ables. Therefore, in an analysis of a Sun-like star, we 
can constrain the fit to all models with this value or use 
a prior based on the marginal posterior probability for 
a m \ as determined from the fit to the Sun. In this case 
we obtain M = 1.04 ± 0.01 M , age = 4.41 ± 0.29 Gyr, 
X = 0.72±0.01, Z = 0.020±0.002, Z s = 0.018±0.001. 
This is an improvement, but still not comparable to the 
results obtained when using the full data set. 

Thanks to the probabilistic method, however, we can 
also easily add new observables as further constraints, 
such as the frequency of maximum power, which can also 
be inferred from a power spectrum analysis and which 
approximately scales for Sun-like stars as 



M/Mq (T eS /T eS>Q 
L/L G 



,3.5 



^max,0 j 



(33) 



with f max ,o = 3120 ± 5 /iHz (Kallinger et al. 2010b). 
Assuming an observed value of kWx.obs and calc ula ting 
^max.mod for each model according to Equation ( 33 ) we 
can then multiply the probability for each model with 



PKax.obslMf,/) 



2-7TOV 



exp 



(^max,obs ^max,mod) 



2a, 2 



(34) 

where a u = y/ 'a 2 (f max ,obs) + 0" 2 (i/ max ,mod)- 

With f m ax.obs = 3120 ± 20 /iHz for our simulated Sun- 
like star, we' then obtain M = 1.02 ± 0.01 M Q , age = 
4.39 ± 0.28 Gyr, X Q = 0.72 ± 0.01, Z = 0.019 ± 0.002, 
Z s = 0.017 ± 0.002. Finally, if we were able to de- 
termine i/jnux b s to about solar precision, the results 
would be M = 1.008 ± 0.006 M , age = 4.39 ± 0.30 Gyr, 
X = 0.71 ±0.01, Z = 0.019±0.002, Z s = 0.017±0.002. 
Therefore, if our observations provide precise additional 
information such as v max , it can easily be implemented 
with our method. It then seems possible to obtain rea- 
sonably accurate results for Sun-like stars, even in the 
absence of low-order modes and without a fixed surface 
effect correction. 

5. CONCLUSIONS 

In this paper, we have derived a new, completely prob- 
abilistic framework for asteroseismic grid fitting. We ex- 
plicitly used marginalization and the formulation of com- 
bined propositions to allow for the quantitative evalua- 
tion of the model grid physics. While computationally 
more intensive than the standard \ 2 evaluation, this ap- 
proach has several benefits in that it 



1. allows for the treatment and analysis of systematic 
errors such as the surface effects, therefore remov- 
ing the need to apply corrections prior to fitting, 

2. easily implements uncertainties in the mode iden- 
tification, 

3. takes into account the fact that grids are discrete 
representations of a continuous parameter space, 
which is especially important for rapidly varying 
bumped modes, 

4. provides a consistent framework to use prior knowl- 
edge about stellar fundamental parameters or to 
evaluate additional observables such as z^ max , and 

5. produces correctly normalized probabilities and 
likelihoods, respectively evidences, which can be 
used to assess the model grid physics and the cali- 
bration of the grids. 

While the above was explicitly derived using the example 
of a static grid, the probabilistic approach would also 
be suited for an adaptive grid approach. The Bayesian 
evidence could be used as a formidable criterion to decide 
whether an adaptive grid needs to be further refined or 
not. 

We also showed how to apply our method to study the 
Sun. The analysis based on our current grid and our 
prior info r matio n matches well the findings of |Houdek| 



& Gough (2011), and in general fits the up-to-date pic- 
ture of the Sun. The age of our best model (measured 
from the pre-main-sequence birth line) is consistent with 
the meteoritic solar age. The solar model arrives on 
the ZAMS approximately 35 ± 5 Myr after appearing on 
the birth line. We found the same best model whether 
non-adiabatic or adiabatic frequencies were used. This 
shows that our method can adequately deal with different 
shapes of surface effects, even when using the same (flex- 
, ible) surface-effect model. One requirement, however, is 
that there exist enough lower-order modes to "anchor" 
the fit. 

To our knowledge, this work is also the first completely 
grid-based asteroseismic analysis of the Sun, using all the 
information provided by the frequencies and prior knowl- 
edge about the solar fundamental parameters, that re- 
sults in the need for initial hydrog en, helium and metal 
mass fractions more consistent with Asplund et al. ( 2009 ) 
than the traditional higher-metallicity models. At least 
for our current grid, these values are required to pro- 
duce a model that "looks" like the Sun, pulsates like the 
Sun, and has the correct solar age. We stress that a for 



mal y fit to th e Sun's oscillation frequencies ( Guenther 



& Brown 20041 or even ta rgeted nonlinear inversio n ot 



the oscillation frequencies ( Marchenkov et aL]|2000 ) will 



not necessarily yield the same model as our approach. 
With x 2 fits it is difficult to provide an unbiased correc- 
tion for surface effects that at the same time does not 
overly weight the deeper penetrating modes. Some of 
the deeper penetrating modes are sensitive to the base of 
the convection zone where the effects of convective over- 
shoot and turbulence, introduced by rotation shears, are 
not included in the standard models. Inversion meth- 
ods, where a standard base model is perturbed to fit the 
oscillations, are also distinct because even though the 
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perturbed model obtained from inversions reveal regions 
of the standard model that are inadequate, e.g., the base 
of the convection zone, the inversion model is not an ac- 
tual standard model in the sense that it is constrained 
and generated by the model physics. 

We know our best-fit model is inaccurate at the surface 
and we suspect it is inaccurate at the base of the convec- 
tion zone (the latter suspicion based on the inadequate 
model physics for this region). Regardless, the model is 
probabilistically the best model from the current model 
grid that matches all the known constraints. We spec- 
ulate that preferring fits that match the oscillation fre- 
quencies at the expense of the other physical constraints 
may be the reason that helioseismologists have been un- 
able reconcile the observed solar p-mode frequencies with 
frequencies derived from standard mode ls based on the 
Asplund mixture and metal abundance (ISerenelli et al. 



[2009||Guzik k, Mussack|2010[ ). We will pursue these mat- 
ters in a future study where we include model grids based 
on the Asplund mixture. 

While the purpose of our analysis of the Sun is to test 
the details of our model physics, our method can also be 
used in general asteroseismic investigations. When ap- 
plying our technique to stars other than the Sun, e.g., re- 
cent asteroseismic targets from the Kepler mission, tight 
prior constraints as in the solar case are generally not 
available. However, the probability formalism can sim- 
ply assign uninformative (e.g., uniform) priors for the 
unknown parameters and still retain all the remaining 
benefits like treatment of missing mode identification and 
of finite grid resolution. 

For current asteroseismology, however, the most im- 
portant feature is the flexible treatment of the surface 
effects that differs from the us ual approach o f empl oy- 
ing the empirical correction by Kjeldsen et al. (2008) to 
the frequencies. Instead of measuring the empirical cor- 
rection for the Sun with the help of a reference model, 
we use a flexible probabilistic model that allows us to 
measure surface effects in any star given our current as- 
teroseismic grids. We do not rely on the validity of the 
solar surface-effect correction and can test new surface- 
effect models that deviate from the solar power-law ap- 
proach. Correctly treating the impact of the surface ef- 
fects on the model probabilities, this also yields correctly 
propagated uncertainties, and therefore a less biased (but 
model-dependent) assessment of the stellar fundamental 
parameters. 

The results presented in the previous section indicate 
that the accuracy of such current asteroseismic analyses 
is still an open question and heavily dependent on the 
number of unaffected, lower-order modes. If there are 
not enough lower-order modes the surface effect will lead 
to systematic errors in the fundamental parameter de- 
termination. However, even in such a case, by looking 
at how the evidence changes as better physics are in- 
cluded in the models, our method can be used to iterate 
toward improved models, hopefully solving the surface- 
effect problem eventually. 
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